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I. INTRODUCTION 

The studies of open quantum systems have a long history There has been renewed 
interest in recent years due to the developments of the conception and possible realization of 
quantum communication and quantum computation . A key concept in the open quan- 
tum system study is the decoherence of a quantum system interacting with environments, 
which plays a very important rule in almost all phenomena in the quantum devices used 
in quantum computation and quantum communication It has been shown that 

the states of an open quantum system will finally relax into a set of "pointer states" in the 
Hilbert space [sl by decoherence, i.e. for a quantum system prepared in a linear superposi- 
tion of its eigenstates, interaction of the system with its environment results in a decay from 
the system's initial pure state, ps(t = 0) = \ipo){ilJo\, to a mixed state, ps(t > 0) = 'YliVipii 
YliPi = 1- To be specific, an arbitrary initial state of the system plus the environment may 
be written as: 

l^(^ = 0)) = ($^C„|n))®|V',), (1) 

n 

where the set \n) stands for the eigenstates of the system and {ipe) is the initial state of the 
environment. This state at time t larger than the decoherence time evolved to a mixed 
state, which may be expanded as: 

I^W) = $^C'm(t)(|m)®|e„)). (2) 

m 

Here, the set of states |m) are the so-called pointer states of the system vAM-, l9|], and le^) 
are the corresponding states of the environment that entangled with |m)'l(l]. A convenient 
way to represent the system interacting with the environment is the reduced density matrix, 
defined as 

where Tr^ means tracing over the environment degrees of freedom. The evolution from ([1]) 
to ([2]) may be rewritten as: 

Ps{0) psit) = J2\Cm{t)\^\m){m\. (3) 

m 

When the time t ^ r^, the non-diagonal elements of the reduced density matrix Ps{t) vanish 
and the diagonal elements achieve their equilibrium values. This effect of decoherence is 
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typical for all known quantum systems that induces an increase of the system's entropy 
and the damping of quantum oscillations with time ll|, Il2[ |. 



A theoretical description of the evolution of the system from tp{0) to il){t) driven 
externally by the environment is generally a very difficult problem. The case that the 
environment is described by Boson fields has been extensively studied in the context of 
Master equation approach, both with Markovian [3] or non-Markovian [l^ approximations. 
Although the master equation scheme can be used for a large number of environments of 

n 

possible types (phonon, photons, etc.) [12], however, the Master equation descri pti on is not 
universally valid for all the models of environment and fragile in some systems 14l |. 



Generally, if the Hamiltonian of the compound system is known, the direct way to solve 
the decoherence problem is to follow the evolution of the compound system over a substantial 
period of time. By setting /i = 1, the time dependent Schrodinger equation is: 

,^ = #^.W. (4) 

Here H is the total Hamiltonian of the system plus the environment. The equation (jl]) 
can be decomposed into a set of first-order ordinary differential equations with the initial 
condition '?/'(0), and the total number of equations is the dimension of the Hilbert space 
of the whole system, which is usually very large. In principle, the set of equations can be 
solved by convenient methods of ordinary differential equations such as Predictor- Corrector 
method or Runge-Kutta method. However, direct solution of the equations will cost too 
much computer resource due to the large number of equations involved. Another scheme 
for propagating equation (jlj) is to expand the evolution operator U{t) = exp{—iHAt) in a 
Taylor series, where At is the time step. 

exp{-iHAt) = 1 - iHAt H . (5) 



It has been stated in Ref. [15[ that a numerical scheme based on this expansion is not 
stable, because it does not conserve the time reversal symmetry of Schrodinger equation. 



Variations of the Tay 



quantum systems 



16 



or series have been proposed and used in calculations of evolution of 



17|. Efficient and stable simulation methods are needed to reduce 



the computation load and to increase the simulation speed. 
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The polynomial expansion method has been used in the calculation of dynamics and/or 
spectral properties of large quantum systems with great success [isl, Q, 19, 20|, Tal-Ezer 
and Kosloff proposed the expansion in terms of the Chebyshev polynomials and tested the 
method with the simple harmonic oscillator and the problem of scattering from a surface, 



high accurate resu 
ventional scheme 



ts were obtained with an efficiency six times higher compared to the con- 



15|,Q. 



Silver and Roder used the Chebyshev polynomial expansion in 



the calculation of density of states of large sparse Hamiltonian matrix 



20| . A fast evolution 



method based on the expans.ou of Chebyshev polynomial for dynamical quantum system. 

was proposed and checked by Loh et al [21]. Dobrovitski et al extended the Chebyshev 
polynomial expansion method in the study of a spin system interacting with a spin-bath , 
obtained the decoherence properties of the system and showed the efficiency and accuracy 
of the method. Since the Chebyshev polynomial is the most frequently used orthogonal 
polynomial in most numerical approximation theories 22|, other kinds of orthogonal poly- 
nomials should also be applicable in the evolution problems. The argument of Chebyshev 
polynomial is bounded to the interval [—1, +1], which is suitable for systems with a bounded 
Hamiltonian, and for systems that only bounded below, a cut off to the energy spectrum is 
inevitable in order to use the method. However, it is well known that some of the orthogonal 
polynomials, like Hermite polynomial and Laguerre polynomial, do not limit their arguments 
to finite intervals. Expansion in terms of these kinds of orthogonal polynomials may have the 
merit in the unbounded systems. In this paper, we will explore the efficiency and accuracy 
of methods based on all these orthogonal polynomials. We constructed methods based on 
the Hermite and Laguerre polynomials and found that the above two mentioned orthogonal 
polynomials do have the required properties. The rest of the paper is organized as follows. 
In Section II, we briefly review the spin-bath model and the difficulty on getting its exact 
solution; In Section HI, three kinds of polynomial scheme will be described for the expansion 
of the evolution operator; In Section IV, we present the results of our numerical simulation; 
Finally, A brief summary is given in the Section V. 

II. HAMILTONIAN 

Two systems are used in this study to test the numerical methods. The flrst is a two 
spin-1/2 system coupled to a spin environment and the second is a particle moving in a 
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double well potential. 



The spin 



reference 



iamiltonian we used in testing our numerical schemes is the one that used in 



2J]. The system consists of two spins-1/2 interacting antiferromagneticly, 
and the system coupled to a bath of non-interacting spins-1/2. The Hamiltonian can be 
written as: 

// = 2JSi-S2 + ^Afc(Si+S2)-Ifc. (6) 

k 

Here Si and S2 are two spins with spin half coupled by the coupling constant J, favoring 
the antiparallel alignment, which constitute the system. The spins Ik, k = 1,2, ■ ■ ■ , N are 
spin half environment spins, interacting with the system by Heisenberg coupling Ak, and 
do not interact with each other. The coupling constant between the two system spins is 
much larger than the couplings to the environment spins, J ^ A^. The couplings A^ are 
uniformly distributed in an interval. Both of the system spins and the environment spins 
can be represented by Pauli matrices. 

The Hilbert space of the whole system is 2^"*"^ dimensional when the environment 
consists of spins. The basis state of the environment can be chosen as the direct product 
of the single states |t) or ||) for each spin J^, here ||) and ||) are eigenstates of the square 
and z components of each spin. For a moderate size of the environment, say, N = 18, 
we have to find an exact solution to about 10^ differential equations. And when A^ is 
increased by one, the number of equations is doubled. For this reason efficient algorithms 
are needed in the studies of the evolution of this kind of problems, especially in the case 
of decoherence where long time simulation is required to reach the pointer state. The 



po 



ynomial expansions based on both Chebyshev polynomial |9 



and Hermite polynomial 



25( 1 are very successful in this case. 



The Hamiltonian for the double well potential is given by: 

H = ^- \uV + \x\ (7) 

where we set m = ^ = 1. This model is very important in the studies of critical phenomena 
and in the standard model of particle physics when the variable x is a scalar field. Here we 
take it to be a simple yet non trivial model to test our numerical method. 



III. POLYNOMIAL SCHEME 



The formal solution of the equation (jlj) is: 

m = e-^^V(O) = ^(t)^(O). (8) 

The evolution operator U(t) is an exponential functional of the Hamiltonian operator H 
which is represented as a matrix in the Hilbert space of ip. The method of polynomial 
expansion is to expand the evolution operator U (t) in terms of the orthogonal polynomials 
of Hamiltonian H. T 
presented in {q] and 



he expansions in Chebyshev polynomial and Hermite polynomial are 



25l | respectively. We will briefly introduce the Chebyshev and Her- 
mite polynomial expansion and give detailed derivation of expansion in terms of Laguerre 
polynomials and check the efficiency of the method numerically. 



A. Chebyshev polynomial 

The Chebyshev expansion of U{t) given by Dobrovitski et al is: 

U{t) = expi-trH) = ^ Ckn{H), (9) 

A;=0 

where r = EQt/2 and H = 2H/Eq, Eq is a scale factor, are the Chebyshev polynomials: 
Tk{x) = cos(A;arccosx). The reason that we change H into H comes from the argument 
domain of Tfc(x), that is x G [—1,1]. For our spins system, H is bounded above and below, 
so that the scale factor Eq can be determined in the following way: 

-Emax = max{tp\H\tp) , 
Emin = mm{tp\H\'ilj) , 

Eo = 2 max( I -Emax, -Emin | ) • 

Using the orthogonal property of Tk, the expansion coefficients Ck of equation ([9]) can be 
calculated as: 

ak Tkexp{-ixT) -Nfc x / ^ 

Ck = — I ^ — dx = ak{-i) Jk{r), 



J-l VT 



X 



where Jk{T) is the Bessel function of the kth order, and ao = 1 when k = and = 2 
when k ^ 1. The series of Chebyshev polynomials of Hamiltonian H can be calculated by 
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the recursion process: 



To{H) = 1, 



B. Hermite polynomial 

In order to obtain the expansion in terms of Hermite polynomials, we start from its 
generating function 26| 

oo ^ 

e"^^^"^ = E|T^^(^)- (10) 

fe=0 

Where Hk{x) denotes the Hermite polynomial of order k. The evolution operator ([8]) can 
be rearranged as 



The second part of the right hand side of equation ( iTTll is identified to be the generating 
function of Hermit polynomial by setting x = XH and s = —it/2X in equation (|TOl) . where A 
is introduced for convenience. From equation (flOi) and (fTTj) we obtain the Hermite expansion 
form of the exponential operator U{t): 



^-iHt ^ g-(t/2A) 



OO / .\U 
2 l-r" 



5^^(t/2A)'=i/,(A^). (12) 

k=0 



The formal solution tp{t) = exp{—iHt)tlj{0) then becomes: 

m = e-(*/^^)^X^^(t/2A)V., (13) 

The Hermite polynomial of H can be obtained by the following recursive algorithm: 

(f>o = V'o, 
01 = 2XHiJo, 
0fc+i = 2XH(j)k - 2k(f)k-i- 
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To discuss the convergence of the expansion, we consider the term when k is large. The 
Hermite polynomial may be replaced by its asymptotical expression 



Hk{x) ^2 2 k^e 2 + 2 ^ cos I v2fc + Ix ^ ) . 



Substitute this into equation f|T2l) and using the Stirling's formula for the factorial, 

A;! ^ exp[A;(lnA; - 1)], A; > 1, 



(14) 



(15) 



the magnitude of the kth term in the expansion of equation f[T2|) for large k is 



{t/2\y 



{t/Xf ^k+i 



2 2 k^e 2" 



knX 



cos 1 V2k + IXH - —J . (16) 



The physically meaningful Hamiltonian should always be bounded below, and for every 
evolution problem, the spectrum of the system has a maximum value determined by the 
initial state, which is in the order of the total energy of the initial state. If we set a 
maximum energy Em, a few times of the total energy, then the states with energy larger 
than this maximum will not enter the calculation, and we have a natural energy cut off of 
the problem, the Em- Then we can replace H in equation f|T6l) with Em to estimate the 
condition of the convergence of the expansion. 



mxy 
k\ 



Hk{XE„ 



fe-i 

< 2"— 



exp 



fc, , k X^E^ 

-\nk-\ \ 

2 2 2 



+ kln 



fc-i j k 
2 2 exp <( -- 

fc-i I k 
2 2 exp <; -- 



In /c — In e + In I — 



-2 



X^E^ 



In 



kX'^ 



X^E^ 



From this expression we see that if 

In 

or the time step t satisfies 

t < 



kX^ 



X^E^ 



> 0, 



A exp 



X^E^ 



(17) 



e - \ 2k 

The kth term is not larger than 2^"^, then the summation is convergent. In the numerical 
calculation given below, we set A = 1/2. 
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C. Laguerre polynomial 



The expansion in terms of Laguerre polynomials can also be derived from its generating 
function |2q |: 

(1 - s)-°-ie^^/(«-i) = ^Ll{x)s\ {\s\ < 1). (18) 

fe=0 

where a distinguishes different types of Laguerre polynomials. By setting s = it/{X + it) 
and X = XH, we get the Laguerre polynomial expansion as: 

The recursion relation of Laguerre polynomials are 

L^ix) = 1 (20) 
^(x) = a + l-x 
{k + l)Ll^,{x) = {2k + a + l-x)Ll{x)-{k + a)L'{,^,{x). 

From the relation we obtain the Laguerre polynomial expansion of Hamiltonian H as: 

00 = m (21) 

(k + 1)0^+1 = (2k + a + l- XHWk -ik + aWk_,. 

Different a gives different choice of the algorithms, the domain of a is in the interval of 
(— l,oo). In the calculation of the spin bath Hamiltonian we use a = — 1/2 and set the 
parameter A = 1 for convenience. For other kinds of Hamiltonian different values of a may 
be used to attain higher efficiency and accuracy. 

The convergency of the expansion of equation (fT9|) is guaranteed by the relationship 



between Laguerre polynomial and Hermite polynomial 



26|: 



L7{^)=^-^H,,{^). (22) 

Substituting equation ([HD, (fT5l) and fl22|) into the expansion term Lj^^^H) and re- 

placing H with Ejn, the total energy of the initial state, we could estimate its asymptotical 
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absolute value by such a procedure: 



it 



1 + it 



t' 



k/2 



2'^2'k'e-'+'^ cos ( V2kTlE^ - 



_\2 _|_ ^2 / 22'te'^(''^'^~-^) 



1 / t 

< 22 



Em 

e 2 



t2 



= exp <|^— A;/2 In 

For large k, and the suitably chosen time step 

'^m + ln2 



l+t^X E^ + ln2 



k 



t < 



exp 



k 



(23) 



the terms approach to zero exponentially. 



It should be noted that the energy cut off Em is only used here for convergence proof. 
In practical calculations, we do not need to specify this cut off and the time step is chosen 
by test and error. 



Comparing to the Chebyshev expansion, the methods of Hermite polynomial and La- 
guerre polynomial have an obvious advantage that no scaling to the Hamiltonian is needed. 
So that these expansions may have wider applications. On the other hand, the recurrence 
relation for both Hermite polynomial and Laguerre polynomial is not numerically absolute 
stable as compared to the recurrence relation of Chebyshev polynomial, which is marginal 



stable 



271 ]. This fact limits the number of terms in the expansion to some value fcmax; the 



effect of numerical instability has little effect for k < k^g^^^ and the effect starts to show up 
beyond this cut off. In the practical calculations the fcmax may be chosen to be 30, and the 
time step is set up accordingly with a specified error tolerance to get convergent results. 
The calculation schemes presented here are very general and are not dependent on the spe- 
cific form of the Hamiltonian, however, the applicability should be tested for each kind of 
Hamiltonian before it can be used in practical simulations. The efficiencies of the three 
kinds of polynomial expansion are almost the same from our numerical calculation, careful 
comparison reveals that for the current models the Laguerre expansion with a = 1/2 is a 
little bit faster than the others. 
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IV. NUMERICAL SIMULATION 



A. Test of the spin model 

The efficiency of_the_Chebyshev expansion over conventional method of calculation has 
already given by 



23|. In this section we checked numerically the efficiency of three 
kinds of polynomial expansions by comparing the performance among the three expansions 
as well as with the predictor- corrector (P-C) and Runge-Kutta (R-K) methods to the spin 
bath Hamiltonian given in section [Tll We calculated two particular variables using the 
Hamiltonian: (i) the 2;-component oscillation of any one of the center spins, i.e. s^, i = 1 
or 2, which demonstrated the decoherence rate of the system; (ii) the time dependence of 
von Neumann's entropy, i.e. = —Trp In p, which characterizes the entanglement degree 
of the state of the system . We use the same parameters as used in jo], Q , the exchange 
strength J = 16.0, are uniformly distributed between and 0.5; The initial condition of 
the system is |'?/'(0)) =|Ti) or written as |10), and the environment is a normalized linear 
superposition of the product states of spins with random coefficients. The time step is 
chosen as At = 0.036, which is determined by the compromise of convergence requirement, 
Tr(p) = 1, and the speed of computation. All of the three schemes are implemented and 
tested, the results are consistent with those given by 



9|, l23|] . We also did the calculation with 
the two widely used ordinary differential equation solver, the predictor- corrector (P-C) and 
Runge-Kutta (R-K) methods. At the request of stability and speed, the time step in these 
two methods is almost 1/10 of that in polynomial scheme. We found that the calculation cost 
of the three polynomial expansion schemes are very close to each other, with the Laguerre 
polynomial expansion is slightly faster, and the results are practically the same. So we only 
give the datum thereinafter obtained by Laguerre polynomial expansion in the following. 

Figure [1] are results of the oscillation of sl(t) and von Neumann's entropy Sy]\f(t) of 
the spin-bath Hamiltonian with parameters given in the figure caption. The results are 
obtained by the Laguerre polynomial expansion method and are consistent with results by 
other methods we tested and those reported in literature j^, 3]. 

The comparison between computation costs of different methods with the same error 
tolerance listed in Table [B From the table we see that: (i) When N is very small, it is hard 
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time time 



(a) Oscillation of Si(t) (b) Evolution of entropy SvN{t) 

FIG. 1: Decoherence of two coupled spins by a spin bath calculated by Laguerre method, the 
parameters are: J = 16, N = 12, the tolerance in obtained this figure is set to be IQ-^. 

TABLE I: Comparison of the R-K method with the polynomial scheme (abbreviated as P-S)for 
the problem of decoherence of spin-bath 



Scheme 


At 


No. 


of hath-spin 


precision 


t 


CPU time 


R-K 


0.0036 




4 


10-6 


9000At 


2 sec 


P-S 


0.036 




4 


10-6 


900At 


2 sec 


R-K 


0.0036 




8 


10-6 


QOOOAt 


406 sec 


P-S 


0.036 




8 


10-6 


900At 


50 sec 


R-K 


0.0036 




10 


10-6 


9000At 


2065 sec 


P-S 


0.036 




10 


10-6 


QOOAt 


242 sec 



to distinguish the calculation speed of the two kinds of numerical computation methods; 
(ii) In general, the speed of polynomial scheme is about 8 times as fast as that of the 
direct solution methods, i.e, the Runge-Kutta (R-K) methods (the corresponding datum 
of predictor-corrector method are almost the same as R-K); (iii) With increasing N, the 
speed advantage becomes more evident. All the data reported here are obtained by a micro 
computer with Intel Pentium M Banias Processor 1400MHz, Memory 256M. 
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B. The double well model with Laguerre polynomial scheme 



A 
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600 700 800 900 1000 



FIG. 2: The time evolution of (x) of three cases: (a), '0(0) = (f)o{x — xq); (b), '0(0) = (f)2{x — xq); 
(c), -0(0) = (f)s{x — xo). All of them are calculated in the condition of A/a; = 0.0013. 



The Laguerre polynomial expansion scheme can easily be extended into the studies of 
continuous quantum systems. As an illustration, we used it in the calculation of the time- 
evolution of a given wave function packet in the double well system. The initial state 
prepared as one of the eigenstates of a harmonic oscillator with unit mass and frequency a;, 
centered at the bottom of the right well, xq = uj/\/4X. That is: 

1/2 



H^{x) is the Hermite polynomial of the mth order. 



Hm{V^{x - Xo)) exp[-a;(a; - Xo)^/2]. 



(24) 



In order to use the Laguerre polynomial expansion scheme in the evaluation of the time 
evolution, we expand the state of the system by a complete basis state. In principle, any 
complete basis can be used in this calculation, however, a better choice of the basis will 
greatly reduce the computation efforts and obtain high accuracy results. In this study we 
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FIG. 3: The period of different \/u with the same initial state V'(O) = 0o(a; — a;0) 

use the eigenstates of a simple harmonic oscillator, (pn{x), n — 0,1, ■ ■ ■ oo, abbreviated as \n) 
as as the expansion basis. The Hamiltonian of the simple harmonic oscillator that defines 



This is not necessarily the optimized basis, however, calculation shows that it is pretty 
good in this problem. 

By introduction of the creation operator a'^ and annihilation operator a, the matrix 
elements of the double-well Hamiltonian can easily be evaluated. The coordinate x and 
momentum p can be represented in terms of the operator and a: 



the basis is 



h = -p + -00 X . 





(26) 
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FIG. 4: The time evolution of Standard deviation of coordinate a = ((x^) — {x)"^)^^"^ of three cases: 
(a), V(0) = (l)o{x — xq); (b), ip{0) = (t>2{x — xq); (c), -0(0) = (t>s{x — xq). Ah of them are calculated 
in the condition of \/uj = 0.0013. 



The action of and a on \n) are: 



a\n 



) = v^|n- 1), 



a^\n) = y/n + l\n + 1), 
1~ 



(27) 
(28) 



h\n) 



uj \ n + 



\n). 



In the a"*" and a representation, the double-well Hamiltonian (I7j) becomes: 

H 



i..[(a+)2 + a2] + A(a+ + a)4. 



(29) 



By using (1271) . the matrix elements of (l29l) can easily be obtained. And the matrix form of 
the Hamitonian can be substituted directly in the Laguerre polynomial expansion scheme 
providing a suitable cut off of the states is specified. In our calculation, we cut off the states 
at n = 49, at which in all cases we studied are already convergent. The initial state ipi^) in 
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the calculation is also expanded in terms of the \n). When m = in (12^ . the expansion is: 



^(0) 



«0 = Xq^I—. 

For other value of m in fl2^ . the coefficients of the expansion can easily be evaluated 
numerically. 

Using the Laguerre polynomial scheme, we calculated the average position (x) and the 
variation cr = ((x^) — (x)^)^^^. Figure [2] plotted the evolution of the average position (x) 
with time. The initial states are the eigenstates of simple harmonic centered at the right 
well of the double well potential. For the state of (j)o{x — Xq), which is located at the Xq 
initially, it oscillates back and forth with time. From figure [2](a) we see clearly the periodic 
motion and the period can easily be identified. The period depends on the value of X/u. 
Smaller X/u corresponding to a deeper well and thus a longer period. Figure [3] plots the 
period as a function of the ratio X/uj, which is decreasing monotonically as expected. For 
the states of higher energies, though the initial state is also localized at the right potential 
well, the average position no longer follows a periodic oscillation between the two wells, 
instead, the particle spends most of the time moving around the center of the potential. 

1 /2 

Figure m are plots of the variation of the position, a = ((x^) — (x)^) , as a function of 
time, which represents the width of the corresponding wave pocket. From the figure, we 
see that for the low energy state (f)o{x — Xq), the width is typically 4 as can be seen in the 
figure, smaller than the total width of the potential at the average energy of (polx — Xq), 
which is about 10, and it looks like a wave pocket bouncing about. The energy of the 
state 4>o{x — Xq) for the parameters chosen is —0.0390, slightly lower than the height of the 
middle peak of the potential. The movement of the center of the particle between the two 
wells is a case of quantum tunneling. In the higher energy cases, the wave pocket spends 
most of the time oscillating around the center of the potential well and there is no well 
defined period can be found. 



A similar problem was studied by Bender et al many years ago 28|]. If we transform the 



X coordinate to q according tog = x + /?/2 and set to = y/KO, the equation is changed 
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into 



(30) 



which is exactly the equation (1) in reference 28|]. We use the same initial conditions as 
used in ^] to calculate (g) by our scheme (Here the number of energy eigenstates is 
truncated to 32, which is sufficient for convergent). The result is given in figure [5l which is 



the same as figure 1 in 



The calculation time for this figure is only about 4 seconds in 



a personal computer of Pentium(R) 4 CPU 2.60GHz, Memory 512M. 




FIG. 5: Time dependence of {q) with (3 = 2.5 



V. SUMMARY 



In summary, we proposed a new polynomial scheme, the Laguerre polynomial expansion 
scheme, and tested its validity and efficiency by means of the spin bath model and a con- 
tinuous double-well model. The obvious merit of this scheme compared to the Chebyshev 
polynomial expansion scheme is that no scaling to Hamiltonian is required, which means 
that a priori knowledge of the lower and upper bounds of the Hamiltonian is not needed. 
On the other hand, the computation efficiency and accuracy of the method are basically 
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the same as the Chebyshev polynomial expansion scheme. 



We have also made use of the Laguerre expansion scheme in other kinds of model 
systems to study the affection of the intra-bath entanglement on the decoherence of the 
center spins. The method is also as efficient and accurate in those models as it was in the 
current spin bath model. The results will be reported in separate publications. 
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Changjiang Scholars and Innovative Research Team in University. 



[1] Ulrich Weiss, Quantum Dissipative Systems (2nd edition), World Scientific 1999. 

[2] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University 

Press, Oxford, 2002). 

[3] M. A. Nielsen and 1. L. Chuang, Quantum Computation and Quantum Information (Cam- 
bridge University Press, Cambridge, England, 2000). 

[4] E. Jogs, H.D. Zeh, et al Decoherence and the Appearance of a Classical World in Quantum 
Theory Springer- Verlag (2003) 

[5] W.H. Zurek, Phys. Rev. D 24, 1516 (1981) 

[6] W.H, Zurek, Rev. Mod. Phys 75, 715 (2003) 

[7] C.W. Gardiner, Quantum Noise Springer- Verlag, Berlin, Heidelberg, New York (1991) 
[8] K. Gottfried, Nature (London) 405, 533 (2000) 
[9] V.V. Dobrovitski, H.A. De Raedt Phys. Rev. E 67, 056702 (2003) 
[10] V. Braginsky and Khalili, Quantum Measurement, Cambridge University Press, London (1992) 
[11] D. Giulini, E. Joos, C. Kiefer, J. Kupsch, I.O. Stamatescu, and H. D. Zeh, Decoherence and 
the Appearance of a Classical World in Quantum Theory Springer- Verlag, Berlin, Heidelberg, 
New York, (1996) 

[12] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P. A. Fisher, A. Garg, and W. Zwerger, Rev. 
Mod. Phys. 59, 1 (1987) 

[13] S. Shresta, C. Anastopoulos, A. Dragulescu, and B.L. Hu, arXiv:quant-ph/0408084 vl 13 Aug 



18 



(2004) 

[14] M. Prasca, Phys. Rev. B 68, 165315 (2003); Phys. Rev. B 71, 073301 (2005) 
[15] H. Tal-Ezer and R. KoslofF, J. Chem. Phys. 81, 3967 (1984). 
[16] A. Askar and S. Cakmak, J. Chem. Phys. 68,2794 (1978). 
[17] D. Kosloff and R. Kosloff, J. Comp. Phys. 52, 35(1983). 

[18] A. Weifie, G. Wellein, A. Alvermann and H. Fehske, Rev. Mod. Phys. 78, 275 (2006). 

[19] R. Kosloff, Annu. Rev. Phys. Chem. 45, 145 (1994). 

[20] R. N. Silver and H. Roder, Phys. Rev. E 56, 4822 (1997). 

[21] Y. L. Loh, S. N. Taraskin, and S. R. Elliott, Phys. Rev. Lett. 84, 2290 (2000); 84, 5028 (2000). 
[22] E.W. Cheney Introduction to Approximation Theory McGraw Hill, New York, (1966) 
[23] V.V. Dobrovitski, H.A. Raedt, M.I. Katsnelson, and B.N. Harmon, Phys. Rev. Lett. 90, 
210401 (2003) 

[24] N.V. Prokofev, P.C.E. Stamp, Rep. Prog. Phys. 63, 669 (2000) 
[25] X.G. Hu, Phys. Rev. E 59, 2471 (1999) 

[26] G. Arfken, Mathematical Methods of Physicists, 3rd ed Academic, New York, (1985) 

[27] W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes, Chapter 

5, Cambridge University Press, Cambridge, (1986). 
[28] CM. Bender and F. Cooper, J.E. O'Dell, L.M. Sinimons.Jr, Phys. Rev. Lett. 55, 901 (1985) 



19 



